Optimal feeding is optimal swimming for all Peclet numbers 
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Cells swimming in viscous fluids create flow fields which influence the transport of relevant nutri- 
ents, and therefore their feeding rate. We propose a modeling approach to the problem of optimal 
feeding at zero Reynolds number. We consider a simplified spherical swimmer deforming its shape 
tangentially in a steady fashion (so-called squirmer). Assuming that the nutrient is a passive scalar 
obeying an advection-diffusion equation, the optimal use of flow fields by the swimmer for feeding 
is determined by maximizing the diffusive flux at the organism surface for a fixed rate of energy 
dissipation in the fluid. The results are obtained through the use of an adjoint-based numerical 
optimization implemented by a Legendre polynomial spectral method. We show that, to within a 
negligible amount, the optimal feeding mechanism consists in putting all the energy expended by 
surface distortion into swimming - so-called treadmill motion - which is also the solution maxi- 
mizing the swimming efficiency. Surprisingly, although the rate of feeding depends strongly on the 
value of the Peclet number, the optimal feeding stroke is shown to be essentially independent of it, 
which is confirmed by asymptotic analysis. Within the context of steady actuation, optimal feeding 
is therefore found to be equivalent to optimal swimming for all Peclet numbers. 



I. INTRODUCTION 

Swimming microorganisms can be found in a variety 
of environments, and encompass a wide range of size and 
locomotion mechanisms (T), Q • For bacteria, motility is 
important to achieve many biological functions, including 
location and migration toward regions rich in nutrients, 
oxygen or light [3L |4, swimming against gravity, or escap- 
ing aggressions 0, Q . Motility is also essential to repro- 
ductive success, in particular for mammals Q • Recently, 
the collective motion of dense swimmer suspensions was 
the focus of a number of studies emphasizing instabilities 
and increased mixing fcj-[l2j]. In order to swim in a vis- 
cous fluid, a microorganism must undergo sequences of 
active and non-time-reversible deformations of its body 
surface [HI, [l4| • This surface deformation sequence will 
be referred to in the following as the stroke, which could 
be either a swimming stroke (leading to a net displace- 
ment of the swimmer center of mass), or non-swimming. 

As such a swimmer performs work against the sur- 
rounding fluid, it creates a flow field and can thus mod- 
ify its immediate environment in an important fashion, 
affecting in particular the transport of nutrients. The 
metabolism of many microorganisms relies on the ab- 
sorption at their surface of various particles or molecules 
which are both diffusing and being advected by the 
swimmer-induced flow. Depending on the organism con- 
sidered, these can range from dissolved gases or low- 
weight molecules, to complex proteins, organic com- 
pounds, small particles, or even sometimes heat. This 
is true from the behavior of small bacteria all the way 
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to large organisms such as the protozoon Paramecium 
which feeds on smaller bacteria, whose typical random 
walk motion is equivalent to a diffusive process at the 
scale of the larger organism fl5l - f]~7| . For simplicity, all 
these cases will be referred to as "nutrients" . 

An interesting transport problem in the dynamics of 
swimming cells concerns the coupling between the flow 
created by the swimmer and the transport of nutrients. 
This coupling can be essential for larger cells or cell 
colonies to achieve feeding rates matching their metabolic 
needs [18[. If n is the diffusivity of the nutrient of inter- 
est, and a the typical scale of the organism, the impact of 
the stroke on feeding is characterized by the value of the 
Peclet number, Pc = Td e f/rdiff, where Td G f is the char- 
acteristic time scale for the shape deformation (stroke) 
and Tdiff = a? / k, is the diffusive time scale around the 
organism. 

At small Peclet number, the concentration gradients 
created by the stroke-induced flow are immediately ho- 
mogenized by diffusion, and therefore shape changes af- 
fect only marginally the instantaneous feeding rate. In 
that case, swimming can still affect feeding indirectly by 
allo wing to access regions of higher nutrient concentra- 
tion [13j | . At large Peclet number, however, the advective 
transport by the flow created by the swimming stroke 
can significantly modify the nutrient concentration field. 
In that case, swimming directly impacts feeding both by 
creating large concentration gradients near the body sur- 
face and by increasing the swimmer ability to scan a large 
volume of fluid 

The purpose of the present paper is to quantify the 
impact of the swimming stroke on the feeding ability of 
an organism and to determine the optimal stroke maxi- 
mizing the nutrient uptake. A priori, the optimal stroke 
should depend on the problem of interest through the 
value of the Peclet number. In Nature, the relevant value 
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of Pc varies by several orders of magnitude, due to the 
large variety of sizes and time scales observed for differ- 
ent microorganisms (from less than 1 /xm for the smallest 
bacteria to several hundred /im for larger eukaryotes) and 
the range of diffusivity coefficients for the nutrients of in- 
terest (in aqueous solutions, kt ~ 10 _ 'ni 2 s -1 for heat, 
K o 2 ~ 3 10~ 9 m 2 s _1 for oxygen and small molecules, and 
K ~ 10 _11 -10 -10 m 2 s" 1 for larger proteins). For a given 
organism, the optimal stroke to maximize, for example, 
heat fluxes might therefore not be the same as the one 
maximizing the absorption of a large protein. 

Performing the swimming stroke has an energetic cost 
for the organism. In this paper we will consider the 
portion of the energy budget which includes the rate 
of working against the fluid, which is instantaneously 
dissipated in the form of heat in the fluid. The or- 
ganism's metabolism imposes a restriction on the max- 
imum energy available for motility, and assuming that 
energy losses other than hydrodynamic can be embedded 
in a fixed metabolic efficiency, optimizing the swimming 
stroke for feeding is a mathematical problem which can 
be formulated as follows: For a given amount of energy 
available to a particular microorganism to create a flow, 
what is the optimal stroke (possibly a non-swimming one) 
that maximizes the nutrient uptake? 

In Stokes flow, both body and fluid inertia are negli- 
gible compared to viscous forces [||. The displacement 
of the microorganism and the hydrodynamic efficiency 
are then entirely determined by the shape change se- 
quence and not by the rate at which this sequence is per- 
formed. Solving for the stroke-induced swimming motion 
and the corresponding velocity field around the organism 
can be tedious for complex geometries, as it generally 
involves the flapping motion of a few or many flexible 
flagella or cilia Most of the available literature fo- 
cuses on two possible modeling approaches. In the first 
one, each flexible appendage or body element is modeled 
individually using slender body theory [l|, [2^ or singu- 
larity methods [21| . The second approach, to which this 
work belongs, considers a simplified geometry for which 
the Stokes equations can be solved exactly. This is the 
case for the classical spherical squirmer model considered 
here [22|, [23| , an envelope model for the dynamics of cil- 
iated microorganisms that has been used previously to 
study h ydr odynamics interactions 0, suspension dy- 
namics [23, [26[ and optimal locomotion |27|. For this 
model, the linearity of Stokes equations can be exploited 
to linearly decompose the stroke in a superposition of 
swimming and non-swimming modes, which can then be 
optimized to maximize the organism displacement for 
a given energetic cost (27l - [30l | . In recent work [27j , we 
showed that the optimal time-periodic swimming strokes, 
i.e. the one leading to the largest swimming speed for 
a given amount of available viscous dissipation, exhibit 
wave patterns reminiscent of the mctachronal waves ob- 
served on the surface of ciliated microorganisms [3l| . 

The effect of swimming on the transport of passive 
scalars has been studied in the past both from Lagrangian 



and Eulcrian points of view. In the Lagrangian approach, 
the capture or drift induced on a given particle by the 
swimming motion of the organism is explicitly solved for 
[l9l l32l [33j . In the Eulerian approach, the organism is 
modeled as being suspended in a continuous concentra- 
tion field of nutrients, and the focus is on the absorption 

flux on the swimmer body ins, mm. 

The feeding of a 

model squirmer was recently addressed for steady and un- 
steady tangential surface motions described by the super- 
position of one swimming and one non-swimming mode 
[34l [H[ . The nutrient uptake was observed to be strongly 
dependent on the value of the Peclet number as well as 
the relative intensity of the non-swimming and swimming 
mode. 

In the current paper we propose to determine the 
optimal feeding stroke for a squirmer, namely the one 
maximizing the uptake of a nutrient by the organism 
for a given hydrodynamic energetic cost. We consider 
the simplest swimmer geometry (a sphere) and focus, 
as our first attempt to solve the problem, on the case 
of a steady stroke where the imposed surface velocity is 
time- independent. Such an assumption is obviously a 
simplification as cilia tips display periodic and unsteady 
displacements. As recently observed [l?}, the optimal un- 
steady stroke for locomotion can in fact be interpreted as 
the periodic regularization of the solution to the steady 
optimal problem. It was also shown ^BJ that for some 
particular limit of large Pe and infinitesimal deformation, 
the average feeding by the unsteady stroke is defined at 
leading order by the result of a modified steady problem. 
Our determination of the optimal steady feeding stroke 
is thus expected to provide important physical insights 
on the relation between swimming and feeding for mi- 
croorganisms. In addition, although results are presented 
here for an idealized organism shape, the optimization 
framework detailed in this paper is applicable to more 
complex geometries and is therefore relevant to a wide 
class of advection-diffusion problems near self-propelled 
organisms. 

In this steady framework, the problem at stake is the 
optimal distribution of the available hydrodynamic en- 
ergy between the different actuation modes of the swim- 
mer, either swimming modes that produce locomotion or 
non-swimming modes that only produce stirring of the 
surrounding fluid. To answer this question, the general 
framework of the steady feeding problem is presented in 
SecHU After considering an organism of arbitrary shape, 
the equations are introduced for the particular case of the 
squirmer and solved numerically for some specific strokes 
using a spectral method, allowing us to gain qualitative 
understanding of the effect of the swimming stroke on 
the concentration field and nutrient uptake. In Sec. IIII| 
we derive an adjoint-based optimization procedure to de- 
termine the optimal stroke for a general swimmer, and 
we apply it to characterize computationally the optimal 
stroke for the squirmer as a function of the Peclet num- 
ber. We show that, to within a negligible quantitative 
difference, optimal feeding is equivalent to optimal swim- 
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ming for all Peclet numbers. Our numerical results are 
compared successfully to predictions of asymptotic anal- 
ysis, at both large and small Peclet numbers. Finally, we 
close by a discussion in Sec. |lVl 



II. NUTRIENT TRANSPORT AROUND A 
SWIMMING MICROORGANISM 

A. Advection-diffusion of a passive scalar near a 
general swimming microorganism 

We consider the transport of a passive scalar field 
around a microorganism which stirs the surrounding fluid 
- and possibly swims as well - by imposing a steady tan- 
gential velocity along its surface, described by u s . The 
surface S and the shape of the organism is therefore as- 
sumed to remain independent of time. Throughout this 
paper, a body-fixed reference frame is considered. The 
Reynolds number, Re = pUa/pi is assumed to be small, 
where U and a are the typical swimming velocity and 
length scale of the swimmer, and p and p, are the density 
and dynamic viscosity of the fluid medium. For Re <C 1, 
fluid and solid inertia can be neglected and the velocity 
field u around the swimmer is solution of the incompress- 
ible Stokes problem 



V f 



u 



-Vp + ^V 2 u = 0, 
V • u = 0, 
u = u s for x£cS, 
— (U + fi x x) for x4co. 



(1) 
(2) 

(3) 
(4) 



In Eq. (|4j) , the translation and rotation velocities, U and 
f2, define the organism swimming motion and are deter- 
mined by imposing the free-swimming conditions of zero 
net hydrodynamic force and torque [36[ 



crndS = 0, 



/ x x (er • n)dS = 0, 
Js 



(5a) 
(5b) 



where er = —pi + /i(Vu + V T u) is the stress tensor in 
the fluid, and n the unit normal vector pointing into 
the fluid. The steady swimming problem in Eqs. (fTj) 
§5§ is linear with respect to u s and its solution for the 
swimming velocities and fluid velocity fields can therefore 
be rewritten formally as 



(U,0) = _Sf u s 



(6) 



where Jz? and C are linear operators depending solely on 
the swimmer geometry. 

The hydrodynamic cost of the swimming motion, V, 
is defined as the rate of work performed by the swimmer 
surface against the fluid, and equal to the energy dissi- 
pation rate by viscous stresses in the entire fluid domain 



V = [ (er : d) dV = - [ u s ■ (er ■ n)dS, 
JVf Js 



(7) 



where d = (Vu + Vu T )/2 is the fluid strain rate tensor. 
In the following, the equations are non-dimensionalized 
using a and ^JVfpa as reference length and velocity 
scales respectively. 

The microorganism is assumed to be suspended in an 
unbounded nutrient solution with concentration C = Coo 
in the far-field. The nutrient is assumed to be totally 
absorbed by diffusion through the swimmer surface, and 
C = is imposed on S. Note that for a real swimmer, 
this assumption is only valid if the nutrient flux at the 
surface is smaller than the cell's metabolic processing rate 
(see Ref. [34| for a discussion of a more realistic boundary 
condition) . 

For convenience, the nutrient concentration is rescaled 
as c = (Coo — C)/Coo. With this rescaling, cC 1 cor- 
responds to near-ambient nutrient concentrations in the 
far-field, while c ~ 1 corresponds to nutrient-depleted 
regions near the organism. The rescaled concentration 
field c(x) is the solution of the steady advection-diffusion 
problem 



Peu • Vc = V^c, 
c = 1 for x£iS, 

c — > for x — > co, 



(8) 
(9) 
(10) 



where u is the velocity field solution of the swimming 
problem in Eqs. U}-©, an d 



(11) 



is the Peclet number defined using the characteristic 
length and velocity scales and the nutrient diffusion con- 
stant, K. 

The flux of nutrient on the swimmer's surface is purely 
diffusive and thus defined in non-dimensional form as 




1 f dC AC 

$ = -— / ^-dS, 



Pe 



dn 



(12) 



where dc/dn = n • Vc and n is the normal unit vector to 
the solid boundary pointing into the fluid domain. When 
u s = 0, the organism (a rigid body) does not create 
any flow field and the energy consumption is V = 0. 
Then Pe = 0, and the nutrient uptake is the solution 
to the purely diffusive problem in Eqs. ([51)- pH|) . with a 
corresponding nutrient flux $o- Rather than the absolute 
nutrient uptake $ resulting from a given stroke, we are 
interested here in its increase relative to the rigid body 
reference case, namely J = 3>/<En = Sh/2, where Sh is 
known as the Sherwood number [34| . 

The problem solved in this paper can be formulated as 
follows. For a given amount of energy available to the 
organism to stir the fluid (measured in a dimensionlcss 
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2n + l 



n(n + 1) 
1 - r 3 



3r 



(17) 



(l-/i^;Wn(r), 
1/1 1 



.(19) 



FIG. 1: Squirmer model and spherical polar coordinates used 
in the paper. On the surface of the swimmer (r = 1), the fluid 
velocity is purely tangential u = Ug(n)ee- In the far-field, 
u ~ —Ue x with U the swimming velocity of the organism. 



The velocity field is easily recovered from ip as 



1 dip 



1 dip 



dp. n/l - \x 2 d\ 



(20) 



fashion by Pe), what is the stroke (i.e. the surface velocity 
field u s ) that maximizes the relative nutrient uptake </? 
Note that non-dimcnsionalizing the problem using the 
energy used by the organism rather than its swimming 
velocity allows for both swimming (U ^ 0) as well as 
non-swimming strokes (U = 0). 



B. The squirmer model 



and the swimming velocity is U = oi\. Using this re- 
lation as well as Eqs. (fl"7|) -(f20" |) . the linear operators C 
and J*f can be expressed in terms of ^-projections on the 
Legendrc polynomials. 

In the squirmer model, the swimming stroke u is thus 
entirely determined by the values of the different mode 
amplitudes a„, or cquivalently the vector a.. Note that 
the non-dimensonalization based on the rate of energy 
dissipation leads to the normalization FiT^ 



The general framework of the previous section is now 
applied to the particular case of a spherical swimmer 
prescribing axisymmetric and steady surface velocities. 
By symmetry, the swimming motion of this so-called 
squirmer is at best a pure translation along a fixed di- 
rection e^, and using spherical polar coordinates with 
respect to this axis centered on the swimmer, all fields 
(velocities, pressure, nutrient concentration) only depend 
on r and fi = cos (9, where is the polar angle with re- 
spect to e x (Figure [T]). By taking a to be the sphere 
radius, the swimmer surface is the unit sphere r = 1, 
and the surface velocity, u s = Ugeg, can be decomposed 
into modes as 1231 



Km) = ^2oi n K n (fi), 



(13) 



with 



KnM= i2n+ ?Z\^ L> M , (14) 
n(n + lj 

where L n (/x) is the n-th Legendrc polynomial. The swim- 
ming stroke is fully characterized by the values of the 
constant coefficients a n (n > 1). The pressure p and 
streamfunction ib can be computed at any point of the 
fluid domain as [H 



p(r,fi) = poo + ^2a n P n (r,n), 



Pn{r,p) = - 



4n 2 -l\ L n (jM) 



n + 1 J r n+1 



(15) 



(16) 



(21) 



with P\ = ai, an d 



2n + 1 a n 
p n = for n > 2. 



(22) 



With this rescaling, all possible strokes correspond to a 
vector (3 on the unit hypersphere (in the remainder of the 
paper, a and (3 will be used equivalently to characterize 
the swimming stroke). Note, that for mathematical con- 
venience, the definition of the Peclet number in Eq. (fTTj) 
was modified to 



Pe 



1 



Va 



K y 127T/1 



(23) 



The particular and so-called "treadmill" squirmer must 
now be pointed out. That swimmer only includes one 
squirming mode (/?„ = S n i) and maximizes the distance 
travelled by the swimmer for a given amount of energy 
[27l [28j . For a general squirmer, the first mode (n = 1) 
entirely defines the swimming velocity, and as such is 
referred to in the following as the swimming mode, as 
opposed to all the other modes (n 7^ 1) which do not 
produce any swimming motion. The second mode (n = 
2) defines the local stress applied by the swimmer on the 
surrounding fluid [1, [H| . 

For a given stroke a, the rescaled nutrient concentra- 
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tion c satisfies the advcction-diffusion problem 



dc d^ n dc 



dr <9/i d/i dr 



c(l,(i) 
c(oo,/i) 



5 fn 2\ dc 



(24b) 
(24c) 



The reference nutrient flux $o corresponds to the case of 
a non-stirring squirmer (i.e. a rigid sphere with a n = 
for all n) for which the solution of Eqs. (|2~4"|) is simply 
c = 1/r. From Eq. (|T2"]) . $ = 47r/Pe, and the relative 
nutrient uptake J takes therefore the simple form 



J 



x dr 



(l,/z)d/z. 



(25) 



C. Numerical computation of the concentration 
field: the Legendre Polynomial Spectral Method 
(LPSM) 

In this section, we outline the numerical method used 
to solve for the advection-diffusion problem, Eqs. (f24)). 
and compute the nutrient uptake for a given stroke ex. 
The method is based on the expansion of the different 
fields using Legendre polynomials in /i and generalizes 
the approach presented in Ref. [lij to the entire stroke 
space. 

The nutrient concentration c(r, fi) is decomposed onto 
Legendre polynomials as 



c(r,n) = E C m (r)L m (n). 



(26) 



m=0 



Substituting Eq. (f2l)| into Eqs. (|2~4"f leads after projection 
on the p-th. Legendre polynomial (p > 0) to a system of 
coupled ODEs in r 



oo oo 



p \ ^ \ ^ / 4 dC m d'i/',, 



rn—0 n—1 



1-mnp dr 

,d 2 C 



(27) 



dC„ 



C P (1) 
CJoo) = 0, 



(28) 
(29) 



where the functions ip n (r) arc defined in Eq. (|19[) and 
^4mnp and B mnp are third order scalar tensors defined in 
Appendix [A] The relative nutrient flux is then obtained 
simply as 



J 



dC 
dr 



(r = l). 



(30) 



In the numerical simulations, the summations in 
Eq. (|27[) arc truncated at a finite number N of squirming 



modes to describe the swimming stroke (1 < n < N), 
and M Legendre polynomial modes are used to describe 
the azimuthal variations of c (0 < m < M — 1). Adapting 
the technique used in Ref. [34[ , the system of ODEs in r 
is discretized on a stretched grid obtained by mapping as 
r = a uniformly-spaced grid of N r points in ^. The 
choice of an exponential stretching allows to cover both 
far-field and near-field concentrations. The function <j) is 
a third-order polynomial in £ such that a fixed fraction 
of the total number of points are contained within the 
expected concentration boundary layer at intermediate 
and high Pe number. The discretized system (|2"7} can 
then be rewritten as 



H C = R, 



(31) 



where C is a N r x M vector containing the values 
of C m (rj) on the different grid points, and H is a 
M x M block-matrix, each block being tridiagonal of size 
N r x N r . The block structure of H is tightly-banded: H 
is diagonal if A = 0, tridiagonal if A = 1, pentadiag- 
onal if N = 2, etc. The contribution to the right-hand 
side R arises from the non-homogeneous boundary con- 
dition on the swimmer surface for the first mode Co(r). 
This large linear system is solved using a direct block- 
Gaussian elimination technique taking advantage of the 
sparse structure of H. 



D. Results 

The Legendre Polynomial Spectral Method (LPSM) 
presented in the previous section is now used to com- 
pute, for different values of Pe, the nutrient concentration 
around a squirmer for simple steady swimming strokes in- 
cluding only the first two squirming modes. In Fig. [21 the 
concentration field is shown for three different values of 
Pe (1, 10 and 100) and three different swimming strokes: 
the pure treadmill (/3„ = <5„i), a combination of modes 1 
and 2, and a pure mode-2 stroke. Note that in the latter 
case, the organism is not swimming. The correspond- 
ing relative nutrient uptake J is given for each case, and 
the streamlines (independent of the value of Pe) are also 
shown for each stroke. 

At low Pe (typically Pe < 1), the concentration dis- 
tribution is close to isotropic, and only a few Legendre 
modes are necessary to compute c(r, fj,) accurately. The 
far- field behavior is reached rapidly, so </>(^ m ax) = 8 is suf- 
ficient with A r ~ 80-100 to achieve errors of at most 0.1- 
0.5% on the nutrient uptake. Note from Fig. [5] that the 
nutrient concentration is not very sensitive to the swim- 
ming stroke, resulting in similar relative nutrient uptake 
J. In that regime, the typical diffusion time is much 
shorter than the advective time, resulting in the homog- 
cnization of the concentration field and a weak front-back 
asymmetry along the swimming direction. 

As the value of Pe is increased, the concentration dis- 
tribution develops a stronger angular asymmetry as a 
nutrient-depleted wake (c close to 1) develops in the 
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J = 1.39 J = 1.1 J =1.0 




FIG. 2: (Color online) Nutrient concentration around the swimmer for Pe = 1, 10 and 100 (from top to bottom) and /32//?i=0, 
5 and oo (from left to right), all the other j3j being taken equal to zero. Far from the swimmer c = 0, while c = 1 at the 
swimmer surface. The dimensionless nutrient flux J is quoted for each case. On the bottom row, the streamlines are displayed 
for each stroke. 



region "behind" the organism. Molecular diffusion is 
not rapid enough to homogenize the sharper advection- 
induced gradients as both processes now act on the same 
time scale. This applies for swimming (/3i 7^ 0) as well 
as non-swimming strokes = 0). 

For Pc 3> 1, a boundary layer develops for the nutrient 
concentration in the region where the flow impinges on 
the swimmer surface. In the regions where the radial flow 
leaves the swimmer surface (wake of the treadmill swim- 
mer or upward direction for the pure stresslet swimmer) a 
nutrient-depleted region forms where molecular diffusion 
processes do not have the time to smooth out the sharp 



concentration gradients induced by the velocity field. Nu- 
merically, more Legendre modes are required (typically 
M ~ 100 for Pe ~ 10 up to M ~ 400 for Pe ~ 400), and 
one needs to extend the r-grid further in the far-field (up 
to </>(£max) ~ 18-20 for the highest values of Pe consid- 
ered) and increase its resolution (up to N r ~ 400 for the 
highest values of Pe considered) . 

Figure [5] shows that for a given stroke the relative nu- 
trient uptake, J, is an increasing function of Pe, empha- 
sizing the systematic benefit of the swimming or stirring 
motion on the feeding process. For a fixed Pe (i.e. con- 
stant energy cost), it also shows that the treadmill swim- 
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mer always performs better than the two other strokes 
considered. The pure treadmill and pure mode-2 strokes 
share the existence of a sharp nutrient-depleted ejection 
zone. However, one notices easily that the gradients at 
the surface of the organism are stronger in the former case 
due to the swimming motion of the organism toward a 
nutrient-rich zone. Swimming appears therefore, in these 
preliminary results, to be a more efficient process than 
simple stirring. The intermediate stroke is characterized 
by a recirculation region that tends to create at high Pe a 
zone of homogenized nutrient concentration in the wake 
of the swimmer, thereby reducing the radial gradients in 
that region as well as the nutrient uptake. 

To confirm this optimality of the treadmill swimmer, 
Fig. [3] shows the value of the nutrient uptake for all 
possible swimming strokes obtained with only the first 
three squirming modes (/3i, /32, P3) (this consists of a two- 
dimensional space because of the constant energy con- 
straint, Eq. (HU). We see in Fig. [3] that, within this 
three-parameter family of flow profiles, the optimal feed- 
ing swimmer is located around /?2 ~ /?3 ~ 0, correspond- 
ing to the treadmill swimmer. 



III. OPTIMAL FEEDING BY A STEADY 
SQUIRMER 

The results of the previous section suggest that the 
treadmill swimmer (/3 n = 8 n i) is a feeding optimum at 
all Peclet numbers. In this section this result is con- 
firmed by considering the formal optimization problem of 
the relative nutrient flux J with respect to the swimming 
stroke a. We start in Sec. lIII Al bv presenting the general 
framework of adjoint-based optimization for a swimmer 
of time-independent arbitrary shape S prescribing tan- 
gential surface velocities u s on its boundaries. Although 
the results presented in the remainder of the paper cor- 
respond to the simplest geometry (a sphere), this frame- 
work can be applied to organisms of arbitrary shapes 
and is of interest for a large variety of advection-diffusion 
problems. We then focus on the particular squirmer con- 
figuration in Sec. IIIIB1 and present our optimization re- 
sults in Sec. IIII CI We show that the optimal feeding 
stroke is essentially the same as the optimal swimming 
one, a result true for all values of the Peclet number. 
These numerical results are also confirmed using asymp- 
totic analysis predictions outlined in Appendices [B] and 

K3 



A. Nutrient uptake gradient for a general swimmer 

To derive the optimal swimmer, the gradient of J with 
respect to the swimming stroke must be mathematically 
determined. This gradient indicates the changes to make 
in the swimming stroke in order to increase J, leading 
to a natural computational implementation of the opti- 
mization search. 



a) Pe = 5 




(6) Pe = 200 




FIG. 3: (Color online) Variations of the relative nutrient flux, 
J, within the (@ 2 , /3 3 )-plane for (a) Pe = 5 and (b) Pe = 200 
is adjusted so that ^2/3j = 1). Nutrient flux isolines are 
also shown in black for clarity and correspond to the values 
indicated in the colorbars. The crosses indicate the position 
of the treadmill swimmer in the (P2, /?3)-plane. 



The gradient is obtained using variational analysis as 
in Ref. [27j. Considering a small variation Su of the 
swimming stroke, and the corresponding change in the 
flow velocity field Su = C ■ Su s (see Eq. 0, the resulting 
change 5*E» is given by 



5<f> 



1 

Pe 



d{Sc) 
dn 



AS, 



(32) 



where n is the outward normal unit vector to the surface 
of the swimmer and, at leading order, 8c is the solution 
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of 



Pcu 



VSc 

5c = 



-\7 2 Sc 
for 



-Pe<5u ■ Vc, 
£ S and x — 



(33a) 
(33b) 



Multiplying Eq. (|33a[) by a test function g and integrat- 
ing over the entire fluid domain Vf, one obtains after 
integration by part that, at leading order, 



(5$ 



/ c{£-5u s ) ■ \7gdV, 
Jv f 



(34) 



provided that the function g satisfies the adjoint equa- 
tion: 

Peu-Vg = -V 2 5 , (35a) 
g = 1 for x G S, (35b) 
g —¥ for oo. (35c) 

Equation (f3"4"| defines the gradient of the absolute nutri- 
ent uptake with respect to the swimming stroke. Since 
$o does not depend on the imposed surface velocity, the 
gradient of the relative nutrient uptake J is obtained sim- 
ilarly. Note that the adjoint field g satisfies the same 
advection-diffusion equation as the original passive scalar 
after replacing Pe by — Pe (or alternatively u s by — u s ), 
so the same analytical or numerical methods can be im- 
plemented to solve for both fields. 



1. At step k, for a given stroke (3^ k \ the LPSM is used 
to solve for the concentration field c and its adjoint 
g. The value of the corresponding nutrient flux 
is also computed from Eq. (|30l) . 



2. From Eqs. fl22J) and fl37|), the gradient J of the 
relative nutrient flux is computed. 

3. At fixed Pe , f3 {k) ■ /3 (fc) = 1 and the gradient tan- 
gential to the unit hypersphere is obtained by pro- 
jection 



V|,J = V P J- ipw.v p J)p 



(38) 



4. VyJ defines the steepest ascent direction on the 
unit hypersphere in /3-spacc and the next iteration 
is carried at a new guess for the optimal f3 
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(fc+i) 



/3 



CO 



sVnJ 



1/3 



(/v) 



(39) 



until convergence is reached to a local maximum 
when it is not possible to find a new guess with 
j(fe+i) > j{k) usm g this procedure, even in the 
limit s — > 0. 



B. Nutrient uptake optimization for a squirmer 

In the particular case of a squirmer, the gradient of the 
relative nutrient uptake J with respect to the swimming 
stroke a is obtained from Eq. ([34]) as 



dJ 



Pe 
"2" 



d^ n dg 3^ n dg 



(36) 

Numerically, both the concentration and adjoint fields 
are determined for a given swimming stroke a using the 
method outlined in Sec. Ill CI The relative nutrient up- 
take J is then obtained from Co(r) as in Eq. ([3H)l . Its 
gradient with respect to a n is computed as 



dJ 

da n 



m— p—0 



A 



rnnp 



dr 



2j> + 1 h 

"2^+1 A c ^ Gmdr 



(37) 



where the functions G m (r) are defined in analogy with 
C m (r) from the adjoint field g{r,fx). All the above inte- 
grals are well defined, taking into account the far-field be- 
havior of c and g and the definitions of A mnp and B mnp . 

In the following, the optimal steady swimming stroke 
for a given energy consumption is determined (i.e. the 
optimal a or (3 at given Pe). Starting from a random 
initial condition on the unit hypersphere, the fol- 
lowing steepest ascent algorithm is applied: 



Results 



1. Optimal squirmer for various Pe numbers 



The preliminary results obtained in Sec. Ill Dl suggest 
that the treadmill swimming stroke corresponds to the 
optimal feeding mechanism at all Peclet numbers. This 
result is confirmed here using the numerical optimiza- 
tion techniques outlined above. An arbitrary stroke is 
characterized by an infinite number of coefficients /3 n ; for 
numerical purpose, this description must be truncated to 
the first N squirming modes, thereby exploring a reduced 
stroke-space. The results of the stroke optimization are 
presented below for the cases N — 3 and N = 8. Com- 
putations performed for larger values of led essentially 
to the same optimal strokes and feeding rates. 

For given values of N and 0.01 < Pe < 300, several 
optimization runs were performed starting with different 
random initial strokes. In each run, a rapid convergence 
was observed toward an optimal stroke, only marginally 
different from the treadmill swimmer (pure mode 1). 

The variation of the optimal feeding rate with the 
Pcclct number, Pe, is shown in Fig. |4^a) and empha- 
sizes the strong gain in feeding rate associated with the 
performance of the swimming and/or stirring motion. As 
,/ = 1 corresponds to the case of a rigid sphere (Pe = 0), 
the quantity plotted on Fig. Ufa), J — 1, is a measure of 
the excess rate of feeding induced by the surface motion. 
Figure BJa) also compares the results of the computa- 
tional optimization procedure for two different values of 
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Appendices [Bl and [Cl (see also Ref. [3^ 
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FIG. 4: (a) Optimal stroke-induced nutrient flux J — 1 and 
(b) Relative difference in nutrient flux, AJ/J, between the 
optimal swimmer and the treadmill swimmer as functions of 
the Peclet number, Pe. Numerical results of the optimiza- 
tion procedure are presented for N = 3 (crosses) and N = 8 
(squares). Several sets of calculations were performed for each 
value of Pe and N. In (a), the solid line corresponds to the 
treadmill swimmer. In (a) and (b), the dashed and dotted 
lines correspond to the asymptotic results for the treadmill 
swimmer at Pe <C 1 and Pe > 1 obtained in Appendices iBl 
and [O 



N with the feeding rate obtained for the treadmill swim- 
mer. The main observation is that although the rate of 
feeding is strongly dependent on the value of the Peclet 
number, the numerical optimal is undistinguishable at 
this scale from that of the treadmill swimmer for all val- 
ues of the Peclet number. The asymptotic scalings for 
the treadmill nutrient uptake Jtreadmiii are obtained in 



^treadmill 
^treadmill 



Pe 

1 + — for Pe < 1, 



2Pc 



for Pe > 1, 



(40) 
(41) 



and show an excellent agreement with the numerical re- 
sults (Fig. g]a). 

The relative difference in nutrient flux, A J/ J, between 
the numerical optimal and that of the treadmill swim- 
mer is shown in Fig. HJb). We see that it is always small 
- below 10~ 3 - across the investigated range of Peclet 
numbers and is maximum around Pe ~ 10. A clear 
power-law scaling can be observed at low Pe for AJ/J; 
for Pe < 1, this power- law behavior is in excellent agree- 
ment with the predictions of the asymptotic analysis (see 
Appendix |Bj : 



AJ 



f 2161 V 
^1034880/ 



Pe 7 w 4.3610" 6 Pe 7 . 



(42) 



As a side note, the computational results above are 
presented only for Pe > 0.05. Below this value, the opti- 
mization algorithm is unable to find optimal strokes per- 
forming better than the treadmill swimmer. This does 
not rule out the existence of a different optimum, but 
indicates that this optimum differs from the treadmill 
swimmer by an amount smaller than the round-off error 
of our computations. 

We now turn to the description of the optimal swim- 
ming stroke. In the /3-space where the stroke is repre- 
sented by a point on the unit hyper-sphere, the differ- 
ence to the treadmill swimmer is measured using the an- 
gle t opt = cos -1 /?i between the optimal stroke (3 and the 
treadmill stroke. The fraction of the stroke energy cost 
not dedicated to the swimming velocity, i.e. the energy 
of the non-swimming modes, is sin 2 t opt and is dir ectly 
related to the swimming efficiency r\ of the stroke [271 ] : 
i op t ~ \J\ — 2rj. The variation of t opt with Pe is shown 
on Fig. [5] For all Pe , the difference between the tread- 
mill and optimal strokes remains small, with t op t < 0.02, 
corresponding to an energy in the non-swimming modes 
accounting to less than 0.04% of the total energy cost. 
However, this small difference depends strongly on Pe 
(see Fig. [5J. At low Pe , t opt scales as Pe 3 , and for Pe < 1 
the numerical results are in excellent agreement with the 
predictions of the asymptotic analysis (see Appendix [Bj : 

2161 o o 
topt ~ ^7440 W °-° 0418 (43) 

For Pe > 1, t op t scales as t opt ~ Pe" 1 / 3 (Fig. [5]). Note 
that the similarity in shape of Figs.^Jb) and [5] is a direct 
result of t opt <C 1, as explained below. 



2. Gradient near the treadmill 

The optimal feeding squirmer is essentially, but not 
exactly, identical to the treadmill swimmer. Therefore, 
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FIG. 5: Dependence on the Peclet number, Pe, of the ori- 
entation angle in /3-space, t opt = cos" 1 /?!, of the optimal 
swimming stroke. As in Fig. [4] results are presented when 
the optimization is performed on N = 3 modes (crosses) and 
N = 8 modes (square). The dashed line corresponds to the 
prediction of the asymptotic analysis at Pe <C 1 obtained in 
Appendix iBl 



its properties and feeding rate are expected to be de- 
termined by the nutrient flux gradient V J in the stroke 
space, evaluated at the treadmill. For t <C 1, the swim- 
ming stroke 



f3 = cosf/3, + sin i /3m 



(44) 



is a perturbation of the stroke from the pure treadmill, 
/3 1 in the direction /3m, such that /3m • f3 1 = (i.e. non- 
swimming stroke). Then, the nutrient flux can be ex- 
panded near f3 1 as 



dJ 



J = Ji +t [ 7777- ) + — 



dp 



t- 




dJ 



(45) 

where derivatives with a 1 subscript are evaluated at the 
treadmill. The nutrient flux is therefore maximum for 
the treadmill if and only if: 



dJ 



and 



d 2 J 



< 



dJ 

W1/1 



(46) 



More generally, the optimal value of t and corresponding 
flux are given at leading order by 



^-opt 

AJ 
J 



(dJ/dfi^- (a 2 J/a/3|[) i 

[(^Wii)i] 2 



2Ji 



(dJ/d/3 1 ) 1 - (d 2 J/d(3 



(47a) 



(47b) 




FIG. 6: Dependence with the Peclet number, Pe, of the nu- 
trient flux gradient dJ/d/3 n with respect to the first four odd 
modes n — 1 (stars), n = 3 (squares), n = 5 (circles) and 
n = 7 (triangles) and evaluated at the treadmill (the even 
mode gradients are equal to zero by symmetry). The power 
law dependence of each component is indicated by a dashed 
line. 



These results emphasize the critical role of the nutrient 
flux gradient dJ/dfij in the localization of the optimal 
feeding stroke with respect to the treadmill. Integrating 
Eq. (|36[) by part, the gradient can be rewritten as 



dJ 

da n 



Pe 



1 j -1 





( dg 


dr 




dc 




dr 


dr J 



dc 
d/i dr. 



(48) 



Using the parity properties in /i of ty n , one easily ob- 
tains that for the treadmill, g(r,fx) — c(r, — fx), and con- 
sequently 



Vp> 1, 



dJ 



dp, 



2p 



= 0. 



(49) 



The variation with Pe of the first four odd (and non- 
trivially zero) gradient components of J at the treadmill 
is plotted in Fig. [6] Clear scalings are identified for Pe <§; 
1 and we obtain by regression [4l| 



0.49 Pc, 



dJ 
My 1 
and more generally 

dJ 



0.0020 Pe 4 



(50) 



dJ 

0(Pe 2p+2 ) farp>l. (51) 

,dp2p+l J j 

These results are in good agreement with the predictions 
of the asymptotic analysis at low Pe (see Appendix [5J 

dJ Pe dJ 2161 4 
d£h~~2' ~ 1034880 6 



0.0021 Pe 4 . (52) 
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Figure [6] also shows that the gradient along the third 
mode always dominates by at least one order of magni- 
tude the gradients in the other directions, for all Pe. The 
difference between the optimal and treadmill strokes is 
therefore expected to be dominated by the third squirm- 
ing mode which is confirmed by the fact that the results 
presented in Section IIII C II are only marginally modified 
between N = 3 and N = 8. 

For large Pe, dJ/d(3 2p +i <C dJ/dj3i, which is consis- 
tent with the asymptotic result that the feeding rate only 
depends on j3\ at leading order (see Appendix |C|) . 



shock-like recovery strokes where material elements on 
the organism surface were brought back to their initial 
position. We conjecture that the same will be true in 
the case of feeding, and that the solution to the optimal 
feeding for periodic surface motion will be a combination 
of the optimal steady (treadmill) with regularization to 
enforce periodicity at a rate allowed by the energetic con- 
straints. Ongoing work in this direction, technically more 
complex as it requires solving for the spatio-temporal 
evolution of both the concentration field and the adjoint 
field, will be reported in the future. 



IV. DISCUSSION 



Acknowledgments 



In this paper, the steady spherical squirmer model 
was used to determine optimal feeding strategies at zero 
Reynolds number. For a nutrient following an advection- 
diffusion equation, we showed computationally and the- 
oretically that, for a fixed amount of energy dissipated 
in the fluid, the optimal feeding mechanism is essentially 
equivalent to the optimal swimming problem, and its so- 
lution maximizes the swimming velocity. 

Perhaps surprisingly, the result that optimal feeding is 
optimal swimming does not depend on the value of the 
Peclet number, which is confirmed by asymptotic anal- 
ysis. At low Peclet, the improvement in feeding rate 
as compared to quiescent fluid environment (pure nu- 
trient diffusion) is, as expected, small and, it increases 
as Pe. This linear scaling arises from the proportional- 
ity between the gain in nutrient uptake and the volume 
swept by the swimming organism, which itself is pro- 
portional to its surface area times its swimming speed. 
In the high Peclet regime, the development of concentra- 
tion boundary layers means that the volume swept by the 
swimming organism decreases, and the relative nutrient 
uptake shows a slower increase with Pe than linear. 

One interesting feature of the optimal feeding (equiva- 
lently, optimal swimming) solution is that it is vorticity 
free. This surface treadmill solution corresponds indeed 
to the only surface distribution of velocity which leads 
to potential flow Stokesian swimming. This result could 
very well be a simple consequence of our emphasis on 
energy cost, as the presence of vorticity always increases 
the rate of energy dissipation [37| ■ 

Note that the occurrence of a Pe-independent optimal 
feeding stroke in our simulations is reminiscent of results 
on optimal tracer mixing obtained for all Pe using flows 
directed from sources to sinks [Hj]. Here, the optimal 
stroke corresponds to the swimmer (a sink) maximizing 
its velocity toward the sources of nutrients in the far-field. 

One of the major assumptions of our modeling ap- 
proach is the restriction of the study to steady surface 
motion. In the case of our work on locomotion optimiza- 
tion [271, we showed that although the treadmill swim- 
mer is itself not physical (due to the non-periodicity of 
the trajectories) , the unsteady optimum was found to be 
a superposition of the treadmill solution with periodic 
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Appendix A: Definition of the A mnp and B mnv 
tensors 



The coefficients A mnp and B mnp used in Section III CI 
are defined in terms of the Lcgcndre polynomials as fol- 
low: 



An 



(2p + l)(2n + l) f 1 



J L m L n L p d/i, (Al) 



D 



(2p + l)(2»+l) f 1 



mnp 



2n(n + 



^ Jjl-l* 2 )L' m L' n L p d». 



They are easily computed using 



.4 



mOp — m p, 



B 



mOp 



0. 



(A2) 



(A3) 



and the following recursive relations for n > 1 



A 

^rnnp — 



2n + 1 



n — 1 
2n-3 



^m,n— 2,p + 



771 + 1 

2m + 1 



B 



mnp 



2m 
2n + 1 
'n(n + 1) 



^ ^-m—l,n—l,p 

(n-2)(n- 1) 



(A4) 



B 



m,n — 2,p 



2n-3 
m(m+l) 

2,m ~\~ 1 in— 1,ti— X,p ■ ti -rn-\-l,n — l.p } 



(A5) 



Appendix B: Asymptotic analyis: optimal feeding 
for Pe < 1 

In this appendix, we focus on the treadmill stroke /3j — 
Sji, and determine the concentration field c, nutrient flux 
J and nutrient flux gradient for Pe C 1 using asymptotic 
analysis. 
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Concentration field around the treadmill 



2. Gradient computation 



For Pe -C 1, c(r, p) is sought in the form of a regular 
perturbation expansion in Pc : 



p=0 



(Bl) 



with Co = 1/r, the rigid sphere (Pe = 0) solution. How- 
ever, this expansion is not uniformly valid over the entire 
fluid domain and one must consider a boundary layer 
at infinity for Pe > [HI, [39| . In the near- field (outer 
solution) , c must satisfy the advection-diffusion equation 



1 



d_ 

Or 



, dc 
dr 



Pe 



or 



J_\ dc 
(B2) 



as well as c = 1 on the swimmer surface. In the bound- 
ary layer ^(p, p) = c(r, p), with p = Per, must instead 
satisfy the boundary-layer equation as 



1 



d_ 
d~ P 



,d^_ 
dp 
\-p? 



d_ 

dp 

dp 



_ Pe 3 



dp 



dtf 
dp 
1 - p 2 d% 



dp 



dp 
(B3) 



as well as ^ — > for p — > oo. 

Both c and ^ are sought as regular perturbation 
series in Pe . Using Matched Asymptotic Expansion 
[ioj . both solutions are computed up to order 0(Pe p ) 
(p = 1,2,3) and integration constants at each order 
are obtained by identifying the two solutions up to 
terms 0(Pe p , Pe p ~ 1 /r, l/r p ) over a matching region 
Pe-p/OH-i) ^r^Pe- 1 . 

The final solution valid up to 0(Pc 4 ) is given by 



c(r, p) 

valid for 1 < r <C Pc -1 and 

f Po 2 
tf(p,p)= Pc+— +Pe 3 



X;Pe^c|(r)L g ( M ), 

p=0 <3=1 



(B4) 



17 + 3£ 
80 8 



x — exp 
P 



{l + p) P 



(B5) 



valid for Pe 1 /^" 1 " 1 ' <c p, where the functions c^(r) are 
defined in Appendix |DJ Using Eq. (|30|). the nutrient flux 
is then obtained as 



J =1 + — 



Pe 
2~ 



13 Pc 
80 



7 Per 
80 



0(Pe 4 ). 



(B6) 



Using the previous expansion and Eq. (|36l) . one can 
compute dJ/d/3 n at the treadmill. Using the front-back 
symmetry of the treadmill velocity field, the asymptotic 
expansion of the adjoint field is obtained as g{r, p) — 
c(r,—p) and &(r,p) = tf(r,—p). Splitting the integral 
in r in Eq. (|36p between inner and outer regions, one 
obtains 



dJ 

da n 

lint 



Pc 



(4 



2 

p e -7/8 



dr 



<J n 



Pe 1 



dp 



' y n 



dg_ 

dp 



dp 



dg 
dr 



dp 



(B7) 

dr 
(B8) 

dp 
(B9) 



where and are projection operators on the Leg- 
endre polynomials 



&n[f](r) = 



2n + l 
n(n + 1) J_ 1 



f(r,p)(l-p 2 )L' n (p)dp, 



^[/](r) = (2n + l) / f(r,p)L n (p)dp 



(BIO) 
(BH) 



For n = 3, the integral in Eq. (JB9J) is at least 0(Pe 4 ). Us- 
ing Eq. (|B4I) and the definition of in Appendix [Ql the 
gradient with respect to the third mode is then computed 

as 

dJ 2161 a , ^ , 
^ = T03i880 Pe + °( Pe >- 

Following a similar approach, the gradient with respect 
to Pi is computed as 



dJ _ Pe 



<3(Pe 2 ). 



(B13) 



Note that a similar but longer approach consists in 
computing the gradients directly from the expansion of 
c for an arbitrary combination of two modes [34] . This 
calculation, omitted here for clarity, also provides the 
second derivative d 2 Jjdfi\ evaluated at the treadmill: 



— = — Pc 2 + OfPe 3 ) 
dp 2 7840 + { 



(B14) 



Using the results of Sec. MI C 21 J has a maximum 
in the (/3i, /^-spacc at /?3 Pt corresponding to a relative 
increase A J/ J of the nutrient flux: 



/3 3 ° Pt 
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Appendix C: Asymptotic analysis: 
at Pe > 1 



optimal feeding 



As shown in Fig. [2j the feeding problem at Pe ^ 1 
is characterized by the formation of a boundary layer in 
the concentration distribution near the squirmer's sur- 
face, whose thickness scales as Pc -1 / 2 due to the bal- 
ance between tangential advection and radial diffusion 
near the swimmer's surface. Generalizing the analysis in 
Ref. [lij to arbitrary strokes, Eq. ([8]) becomes at leading 
order in 1/vPe , 



(CI) 



with R — a/Pc (r — 1) and £ the axial component of the 
tangential surface velocity: 



CM = y/l-H 2 u$(li). 



(C2) 



A self-similar solution is sought for the previous equa- 
tion in terms of the variable rj = R/g(p), where <?(//) 
represents the boundary layer thickness. Equation (|C1[) 
then becomes 

eft" c f)c 

W = J1 dv tetogWW) + .9(m) 2 C'(^)} ■ (C3) 



Provided that 

C(^)5W(M)+5(M) 2 C'(^) = -2, 



(C4) 



a self-similar solution compatible with the boundary con- 
ditions Eq. ([9|)- ([T0|) exists, given by 



c(R,n) = —= 



dry. 



(C5) 



Equation (|C4j) can be solved for g{fj) with the addi- 
tional constraint that the boundary layer thickness g(fi) 
is finite at /i = 1 



SO) 




(C6) 



If the surface velocity Ug is positive everywhere (i.e. there 
are no recirculation regions), then g(fJ-) remains finite for 
all /i > —1, and the boundary layer solution above is valid 



over the entire surface of the squirmer. This condition is 
satisfied by the treadmill swimmer, and in some vicinity 
of it. At /i = — 1, g(p) = oo and the boundary layer 
separates into the wake observed on Fig. [5J 

The nutrient flux at the surface of the squirmer is then 
obtained from Eqs. (j2"5)) and (|C5[) as 



J 



d/i 



n J-i 9(jj) 

and can be computed exactly from Eq. (|C6j) as 



(C7) 



J = 




C(/*)d/* + 0(l) = 



2/JxPe 



0{l). (C8) 



One observes that at leading order J depends on /3i only. 
Appendix D: Definition of the functions c|(r) 



The functions c2(r) in Eq. (|B4[) are given by 
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